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Abstract 

Direct  numerical  simulations  (DNS)  of  spanwise-rotating  turbulent  channel  flow  as  well  as  the 
neutral  and  unstable  turbulent  Ekman  layer  were  conducted.  These  DNS  results  were  used  to 
evaluate  various  turbulence  and  heat  transfer  models  for  the  Reynolds  stresses,  turbulent  heat 
fluxes  and  higher-order  moments  of  velocity  and  temperature.  Explicit  Algebraic  Reynolds  Stress 
Models  (EARSM)  obtained  the  Reynolds  stress  distributions  in  best  agreement  with  DNS  data  for 
rotational  flows  and  turbulent  heat  flux  distributions  obtained  from  two  explicit  algebraic  heat  flux 
models  consistently  displayed  increasing  disagreement  with  DNS  data  with  increasing  rotation  rate. 
DNS  results  were  also  used  to  determine  the  proper  computational  box  size  for  a  minimal  flow  unit 
(MFU)  at  i?Ob  =  0.5,  spanwise  arrays  of  Taylor-Gortler  vortices  in  the  highly  turbulent  pressure 
region  were  examined  and  complete  realization  of  the  vortices  was  demonstrated  to  be  necessary 
for  accurate  MFU  turbulence  statistics  requiring  a  minimum  spanwise  domain  length  Lz  —  n.  For 
the  neutrally  stratified  Ekman  layer,  the  higher-order  moments  of  velocity  were  examined  and 
the  accuracy  of  a  kurtosis  model  was  assessed.  For  the  unstable  Ekman  layer,  the  analysis  of 
higher-order  moments  was  extended  to  temperature- velocity  correlations.  Model  coefficients  were 
optimised  using  DNS  data  and  it  was  shown  that  the  optimised  models  accurately  captured  the 
distributions  of  all  fourth-order  moments.  These  flow  fields  represent  complex  turbulence  which 
will  be  subject  to  flow  and  heat  transfer  control  by  phonons  at  a  later  stage  of  our  work.  Research 
aimed  at  the  control  of  fully  turbulent  channel  flow  using  direct  numerical  simulation  (DNS)  was 
also  conducted.  The  reduction  of  the  kinetic  energy  of  large  amplitude  perturbations  in  channel 
flow  was  investigated  using  passive  phononic  (periodic)  structures.  These  studies,  and  the  results 
obtained,  lay  the  foundation  for  extending  the  phononic  subsurface  passive  control  methodology 
to  turbulent  drag  in  channel  flows.  Results  from  this  seed  grant  have  appeared  in  three  journal 
articles  and  have  been  accepted  for  presentation  at  national  conferences. 
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I.  RESEARCH 

Several  different  research  endeavors  were  conducted,  all  involving  the  integration  of  the 
time-dependent  Navier-Stokes  equations  using  direct  numerical  simulation  (DNS). 


A.  Spanwise-Rotating  Turbulent  Channel  Flow  and  Modeling 

Turbulent  channel  flow  subject  to  rotation  in  the  spanwise  direction  is  characterized  by 
reduced  turbulence  levels  near  one  wall  and  elevated  turbulence  levels  near  the  opposite  wall; 
these  regions  are  known  as  the  suction  and  pressure  sides,  respectively1.  Subsequently,  the 
symmetric  profiles  of  mean  velocity  and  Reynolds  stress  distributions  in  the  non-rotating 
channel  become  asymmetric  with  respect  to  the  channel  centerline.  The  effects  of  spanwise 
rotation  on  momentum  transport  in  turbulent  channel  flow  has  been  well  documented  in 
previous  direct  numerical  simulation  (DNS)  studies  by  Grundestam  et  al.1  and  Wu  and 
Kasagi2.  Rotation- induced  body  forces  (Coriolis,  centrifugal)  generate  a  secondary  cross 
flow,  and  consequently  a  particular  type  of  complex  turbulent  flow  regime  is  developed  with 
more  than  one  mean  flow  gradient.  The  analyses  of  such  complexities  on  the  structure 
and  parameterization  of  turbulence  have  relevance  to  engineering  applications  such  as  gas 
turbine  blade  and  rotating  turbomachinery  design,  especially  with  regards  to  surface  heat 
transfer  and  skin  friction  within  the  internal  cooling  passages3,4. 

For  the  research  conducted  and  published  in  Hsieh,  Biringen,  and  Kucala5,  DNS  was 
employed  to  assess  four  RANS  models  proposed  by  (a)  Reif  et  al.6  (PRDO),  (b)  Speziale 
and  Gatski7  (SG),  (c)  Girimaji8  (GI)  and  (d)  Grundestam  et  al.9  (GWJ).  Two  algebraic  heat 
flux  models  proposed  by  (e)  Younis  et  al.10  (YWL)  and  (f)  Abe  and  Suga11  (SA)  were  also 
evaluated.  In  addition,  the  pressure-strain  functions  proposed  in  Speziale  and  Gatski7  and 
Girimaji8  were  investigated  for  their  influence  on  the  modeled  Reynolds  stress  distributions. 

For  spanwise-rotating  turbulent  channel  flow,  the  Reynolds  stress  distributions  produced 
from  a  linear  and  nonlinear  eddy  viscosity  model  were  compared  to  demonstrate  improved 
accuracy  for  the  nonlinear  model  over  the  linear  model.  In  a  comparison  of  four  nonlinear 
eddy  viscosity  models  with  DNS  data,  EARSM  were  the  most  compatible  with  DNS  results 
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in  modeling  the  Reynolds  stresses  for  turbulent  channel  flow  subject  to  spanwise  rotation. 
The  Speziale-Gatski  model  was  shown  to  be  the  most  compatible  for  zero  and  low  rotation 
numbers  but  displayed  significant  deviations  near  the  pressure  wall  at  high  rotation  numbers. 
As  shown  in  figure  1,  the  Grundestam- Wallin- Johansson  model  showed  the  best  agreement 
with  the  DNS  data  at  high  rotation  numbers.  Heat  flux  models  were  in  good  agreement 
with  DNS  data  for  the  no-rotation  case  but  with  system  rotation,  the  models  deviated  from 
the  DNS,  increasing  at  higher  rotation  rates. 


FIG.  I:  Modeled  Reynolds  stress  profiles  for  Ro^  =  0.9.  (a)  upjj  ;  (b)  u2u2  ;  (c)  u3u3  ; 
(d)  u{u2  •  Black:  DNS;  blue:  PRDO;  green:  SG;  red:  GI;  magenta:  GWJ. 


The  pressure-strain  models  of  two  EARSM  (Girimaji,  SG)  were  shown  to  have  significant 
disagreements  with  the  DNS  data  in  the  near- wall  regions  and  the  pressure-temperature- 
gradient  models  of  two  EAHFM  (YWL,  SA)  demonstrated  inaccurate  characterization  of 
the  suction  region  with  system  rotation.  The  errors  in  the  modeled  contributions  from 
these  terms  resulted  in  degeneration  of  the  predictive  capabilities  of  their  respective  closure 
models.  These  errors  contributed  to  inaccurate  Reynolds  stress  amplitudes  in  the  near- wall 
regions  and  an  inaccurate  modeled  distribution  shape  for  wall-normal  turbulent  heat  flux 
in  EARSM  and  EAHFM,  respectively.  Present  results  indicate  correct  characterization  of 
pressure  fluctuations  is  a  crucial  factor  in  both  EARSM  and  EAHFM  design  for  spanwise- 
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FIG.  2:  Three-dimensional  contours  of  time- averaged  v  and  w  velocity  for  Ro\,  =  0.5.  Dark 
and  light  contours  denote  clockwise  and  counter-clockwise  motion,  respectively. 

rotating  turbulent  channel  flow. 


B.  Spanwise-Rotating  Turbulent  Channel  Flow  and  the  Minimal  Flow  Unit 

For  spanwise-rotating  turbulent  channel  flow,  the  streaky  and  vortical  structures  associ¬ 
ated  with  the  turbulence  sustenance  cycle12,13  persist  in  the  pressure  region  and  the  gener¬ 
ation  of  additional  turbulence  structures  was  observed  in  the  rotational  turbulence  studies 
of  Kristofferson  and  Andersson14  and  Grundestam  et  al.1 .  In  figure  2,  one  example  of  the 
rotation-induced  flow  instabilities,  known  as  roll  cells,  are  shown  to  circulate  flow  throughout 
the  pressure  region  of  the  chnnel. 

In  order  to  determine  the  dimensions  of  a  minimal  flow  unit,  it  is  imperative  to  examine 
the  contributions  of  these  rotation-induced  structures  to  turbulence.  The  concept  of  the  min¬ 
imal  flow  unit  (MFU)  model  is  based  on  the  determination  of  the  smallest  computational  box 
size  that  will  produce  acceptably  accurate  turbulence  statistics  at  minimal  computational 
cost.  As  the  dependence  of  turbulence  production  on  the  interactions  of  various  turbulence 
structures  has  been  well  documented  in  literature13,15,  MFU  design  distinguishes  a  basic 
set  of  structures  necessary  to  sustain  turbulence  and  constructs  a  shorter  domain  based  on 
this  array  of  structures.  A  model’s  success  is  determined  by  its  ability  to  accurately  predict 
essential  turbulence  statistical  quantities  at  a  significantly  reduced  computational  cost  com- 
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FIG.  3:  Three-dimensional  contours  of  time-averaged  v  and  w  velocity  for  MFU  simulation 
case  at  R(>\,  =  0.5.  Dark  and  light  contours  denote  clockwise  and  counter-clockwise  motion, 

respectively. 

pared  to  full  direct  numerical  simulation  (DNS)  on  a  well-resolved  computational  domain. 
Such  MFU  models  are  necessary  for  computational  fluid  dynamics  research  requiring  large 
amounts  of  simulations  such  as  parametric  studies  for  flow  control16,17.  For  the  research 
conducted  and  published  in  Hsieh  and  Biringen18,  the  DNS  database  for  turbulent  channel 
flow,  subject  to  varying  rotation  and  Reynolds  numbers,  was  used  to  assess  the  accuracy  of 
the  MFU  model  for  predicting  low  and  high-order  moments  of  turbulent  fluctuations. 

For  the  design  of  a  minimal  flow  unit  for  rotational  turbulence,  a  baseline  MFU  model 
with  spanwise  domain  length  of  Lz  =  ir8  was  selected  to  accomodate  a  single  full  pair  of 
Taylor-Gortler  vortices  as  shown  in  figure  3.  A  box  minimization  study  with  reduced  span- 
wise  domain  lengths  down  to  Lz  =  0.187T<5,  the  MFU  length  for  the  non-rotating  turbulent 
channel  flow,  was  conducted.  Observed  discrepancies  in  the  mean  velocity  distributions 
demonstrated  that  MFU  accuracy  did  not  depend  on  sublayer  streak  distance  as  for  the 
non-rotational  channel  and  a  significantly  larger  minimum  spanwise  length  Lz  =  nd  was 
required  for  accurate  turbulent  statistics,  corresponding  to  the  minimum  length  for  proper 
realization  of  one  full  pair  of  Taylor-Gortler  vortices.  If  these  vortices  were  inaccurately 
represented  from  further  truncation  of  the  spanwise  domain  length,  turbulent  fluctuations 
were  inaccurate  or  an  incorrect  mean  velocity  gradient  was  produced  in  the  pressure  region. 
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For  a  higher-Reynolds  number,  the  MFU  model  demonstrated  decreased  accuracy  com¬ 
pared  to  the  low-Reynolds  simulation.  Hence  for  large  Reynolds  numbers,  MFU  models 
may  require  a  significantly  larger  domain  box  to  accurately  approximate  turbulence  statis¬ 
tics  and  alternative  factors  for  MFU  design  require  consideration,  such  as  Reynolds  number 
effects  on  sublayer  streak  length  and  turbulence  structures  in  the  suction  region.  To  test 
the  limitations  of  the  MFU  model,  higher-order  statistics  from  the  baseline  MFU  model 
were  compared  to  those  from  the  full  simulations.  The  model  produced  accurate  distri¬ 
butions  of  skewness  and  kurtosis  for  a  non-rotating  channel  but  was  unable  to  maintain 
this  accuracy  with  rotation  in  the  suction  region.  The  MFU  model  accurately  captured 
higher-order  statistics  in  the  pressure  region  due  to  the  successful  realization  of  roll  cells  but 
could  not  properly  capture  the  re-laminarized  suction  region  which  contained  intermittent 
high-amplitude  velocity  fluctuations,  a  consequence  of  the  turbulent  spots  structures.  These 
findings  indicated  that  when  the  MFU  model  was  extended  beyond  its  intended  function 
of  general  turbulence  quantities  (mean  velocity,  Reynolds  stresses)  to  higher-order  statis¬ 
tics,  the  model  continued  to  perform  well  in  regions  of  high  turbulence  due  to  its  ability  to 
capture  the  coherent  structures  which  contribute  to  turbulence  production.  However,  the 
model  failed  in  regions  with  different  physical  dynamics  such  as  the  low-turbulence  suction 
region. 


C.  Ekman  Layer  Flow  and  Modeling 

The  Ekman  layer19  is  a  boundary  layer  formed  by  pressure  gradients  in  a  rotating 
system20.  With  a  heated  surface  (convective  boundary  layer)  and  capped  inversion,  the  Ek¬ 
man  layer  is  often  used  to  model  the  complex  dynamics  of  the  atmospheric  boundary  layer 
(ABL)  such  as  buoyant  forcing  and  effects  of  Coriolis  forces  due  to  the  Earth’s  rotation21,22. 
Similar  to  the  turbulent  channel  problem,  turbulence  closure  models  that  consider  time- 
averaged  equations  with  phenomenological  closure  approximations  are  highly  desirable  if 
they  can  accurately  parameterise  turbulent  transport. 

For  the  research  conducted  and  published  in  Waggy,  Hsieh,  and  Biringen23,  the  DNS 
database  under  both  neutral  and  unstable  conditions  was  utilized  to  assess  closure  models 
used  for  predicting  high-order  moments  of  turbulence  and  temperature  fluctuations  as  a 
function  of  lower-order  correlations.  The  higher-order  moments  of  skewness  and  kurtosis 
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in  the  turbulent  Ekman  layer  were  examined;  an  assessment  of  a  model  proposed  by  Mole 
and  Clarke24  for  the  kurtosis  was  also  provided.  For  the  unstable  Ekman  layer,  analysis  of 
higher-order  moments  was  extended  to  temperature- velocity  correlations,  and  two  closure 
models  by  Zilitinkevich  et  al.25  and  Gryanik  and  Hartmann26  were  evaluated. 

Evaluation  of  the  DNS  results  with  previous  similar  studies27,28  showed  strong  agreement, 
lending  validity  to  the  simulation  results  used  for  turbulence  modeling.  The  higher-order 
moments  of  skewness  and  kurtosis  were  introduced  in  the  case  of  a  neutrally  stratified 
Ekman  boundary  layer  and  a  model  proposed  by  Mole  and  Clarke24  for  the  kurtosis  was 
evaluated.  An  analysis  of  two  separate  sets  of  coefficients  proposed  for  the  model  by  Tampieri 
et  al. 29  and  Alberghi  et  al. 30  demonstrated  that  the  coefficients  proposed  by  Alberghi  et  al. 30 
generally  had  better  agreement  with  the  DNS  data.  In  the  case  of  an  unstratffied  Ekman 
boundary  layer,  two  closure  models  by  Zilitinkevich  et  al.25  and  Gryanik  and  Hartmann26 
approximated  third  and  fourth-order  moments  of  velocity  and  temperature  as  a  function  of 
lower  order  moments  as  well  as  the  skewness  and  kurtosis.  A  parametric  study  using  the 
explained  variance  (a2)  was  conducted  to  assess  the  accuracy  of  the  coefficients  proposed 
by  each  of  these  models  and  propose  a  new  set  of  coefficients  that  would  maximise  a2.  The 
computed  model  coefficients  that  maximised  a2  did  well  in  accurately  capturing  the  trend 
of  all  fourth-order  moments  as  well  as  some  third-order  moments. 


D.  Flow  Control:  Channel  Flow 

Flow  control  in  regards  to  jet  turbines  has  been  a  subject  of  great  interest  in  recent  years. 
Turbulence  is  an  impediment  to  effective  turbine  design  due  to  its  flow  regime  consisting 
of  chaotic  property  changes  and  instabilities,  leading  to  undesirable  results  such  as  higher 
drag  and  energy  losses31.  Flow  control  is  also  a  major  topic  of  research  for  friction  drag 
reduction  which  too  is  directly  connected  to  instabilities32.  Flow  control  focuses  on  the 
reduction  of  these  instabilities  using  active  (non-zero  energy  cost)  or  passive  (zero  energy 
cost)  methods.  This  research  examined  passive  flow  control  using  phononic  crystals  placed 
underneath  the  surface.  A  phononic  crystal  is  a  periodic  material  formed  by  the  repeated 
spatial  arrangement  of  a  unit  cell33,34.  The  unit  cell  exhibits  a  band  structure  that  relates 
the  frequency  to  the  wavenumber  of  elastic  waves  traveling  across  the  phononic  crystal  as 
a  whole.  Among  the  features  of  a  band  structure  that  is  widely  used  in  phononic-crystal 
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Design  A 

(PRSA) 


FIG.  4:  Wavenumber,  amplitude,  phase  and  performance  metric  curves  vs.  frequency  for  a 

phononic  subsurface  design  labeled  Design  A. 


applications  is  the  presence  of  band  gaps  (frequency  ranges  were  waves  are  prohibited  from 
propagation35.  In  the  context  of  flow  control,  our  previous  research  has  demonstrated  that 
within  a  stop  band,  not  only  the  waves  are  prohibited  from  propagation,  but  their  phase 
changes  in  a  robust  manner,  i.e.,  independent  of  the  boundary  conditions.  When  combining 
this  effect  to  the  tuning  of  the  finite  phononic  crystal  resonance  response,  we  have  shown 
that  it  is  possible  to  use  this  mechanism  for  passive  stabilization  of  a  laminar  flow  exhibiting 
a  growing  Tollmien-Schlichting  (TS)  instability17. 

As  shown  in  figure  4,  the  effect  of  a  particular  phononic  subsurface  design  over  a  large 
frequency  range  varies  significantly  depending  on  the  layering  of  individual  materials  within 
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FIG.  5:  Streamwise  spatial  distributions  of  the  kinetic  energy  of  the  disturbance  field 
within  the  bottom  half  of  the  channel;  the  dotted  lines  represents  the  base  model  (i.e., 
with  no  control)  and  the  continuous  lines  represent  the  model  with  the  phononic  crystal. 
The  plotted  quantity  represents  the  spatial  intensity  of  the  flow  instability. 


the  unit  cell.  Despite  the  material  system  as  a  whole  being  composed  of  only  two  constituent 
materials  (ABS  polymer  and  aluminum),  the  number  of  layers  as  well  as  the  volume  fraction 
of  each  constituent  material  may  be  modified  freely.  Hence  optimization  studies  can  be  per¬ 
formed  to  tuning  these  phononic  subsurfaces  to  suppress  the  turbulence  mechanisms.  Our 
research  aims  to  produce  an  optimized  phononic  subsurface  structure  which  significantly 
reduces  both  turbulent  kinetic  energy  and  drag  in  turbulent  channel  flow. 

Prior  to  investigating  the  turbulence  problem,  it  is  important  to  examine  the  effects 
of  large  amplitude  disturbances  since  these  bring  rise  to  nonlinear  effects.  To  examine  the 
ability  of  a  phononic  subsurface  to  counter  large- amplitude  instabilities,  we  conducted  cou¬ 
pled  fluid-structure  direct  numerical  simulations  where  we  introduced  an  initial  excitation 
wave  incorporated  as  a  spatially  evolving  disturbance  in  a  fully-developed  plane  Poiseuille 
(channel)  flow  driven  by  a  mean  pressure  gradient.  The  channel  is  formed  by  parallel  walls 
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at  the  bottom  and  at  the  top  with  periodic  boundary  conditions  applied  in  the  spanwise 
^-direction,  and  a  buffer  (sponge)  layer  is  used  to  model  the  outflow.  In  the  fluid  part  of 
the  coupled  model,  the  base  flow  is  an  exact  solution  of  the  Navier-Stokes  equation.  We 
superimposed  an  unstable  spatial  solution  (eigenfunction)  of  the  Orr-Sommerfeld  equation 
at  the  inflow  boundary  of  the  channel  to  excite  the  parabolic  base  velocity,  faithfully  mod¬ 
eling  the  conditions  in  typical  laboratory  experiments.  We  consider  water  as  the  working 
fluid,  a  Reynolds  number  Re  =  7500,  and  a  non-dimensional  unstable  frequency  ujr  =  0.25 
(i.e.  1690Hz).  In  the  solid  domain,  the  Newmark  scheme  was  used  for  the  time  integra¬ 
tion.  The  unit  cell  of  the  phononic  crystal  configuration  considered  was  again  composed  of 
aluminum  and  ABS  polymer.  The  density,  ps  and  Youngs  modulus,  Es,  for  each  of  these 
two  constituent  materials  are:  pai  =  2700  Kg/m3,  Pabs  =  1040  Kg/m3,  Eai  =  68.8  GPa, 
Eabs  =  2.4  GPa.  Both  the  structure  and  the  fluid  equations  are  inverted  separately  and  a 
conventional  serial  staggered  approach  is  used  to  couple  the  interface  between  the  two. 

The  streamwise  evolution  of  selected  modal  contributions  of  perturbation  kinetic  energy 
(KE)  is  shown  in  Fig.  5(a).  Compared  to  the  reference  all-rigid- wall  case  (dashed  lines), 
it  is  observed  that  the  perturbation  KE  decreases  across  the  length  of  the  phononic  crystal 
interface  for  the  primary  mode  (blue).  However,  this  effect  is  reversed  for  the  three  other 
modes  presented  and  may  be  explained  by  nonlinear  interactions.  This  reversed  effect  is 
somewhat  negligible,  however,  owing  to  the  difference  in  the  orders  of  magnitude  of  the 
modal  energies.  By  only  stabilizing  the  primary  mode,  we  show  a  reduction  of  maximum 
8%  in  the  total  perturbation  kinetic  energy  (summed  over  all  modes)  as  shown  in  Fig. 
5(b).  These  results  demonstrate  the  ability  of  the  phononic  crystal  to  stabilize  the  three- 
dimensional  flow  field  [the  type  exhibited  by  flows  undergoing  secondary  (K-type)  transition] 
via  frequency-dependent  wave  interferences,  even  in  the  presence  of  large-amplitude,  non¬ 
linear  disturbances. 


II.  COMPUTATIONAL  METHOD 

Direct  simulations  are  important  to  the  scientific  community  because  of  the  detail  they 
offer  into  the  dynamics  of  a  flow.  Near  boundaries  a  fine  mesh  is  required  to  capture 
the  smallest  scales  of  turbulent  motion.  In  many  instances,  DNS  allows  monitoring  of 
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FIG.  6:  Diagram  of  turbulent  channel  flow  with  rotation  in  the  spanwise  direction. 


flow  variables  closer  to  the  wall  than  experimentation.  DNS  is  the  ideal  tool  to  use  in 
problems  dealing  with  flow  control,  where  very  precise  calculation  of  the  near-wall  dynamics 
is  critical.  Other  methods,  such  as  Reynolds  Averaged  Navier-Stokes  (RANS)  and  large- 
eddy  simulations  (LES),  while  able  to  capture  the  large  scales  in  the  problem,  are  unable  to 
resolve  smaller  scales  of  turbulence.  The  use  of  DNS  bypasses  this  problem,  but  at  a  much 
higher  cost  of  computational  resources. 


A.  Numerical  Algorithm  and  Parallel  Architecture 

For  the  problems  described  above,  we  incorporate  a  semi-implicit  finite  difference  method 
to  solve  the  incompressible  Navier-Stokes  (N-S)  equations.  The  nonlinear  advective  terms 
are  solved  using  a  second-order  time  and  fourth-order  spatial  variant  of  the  Adams-Bashforth 
explicit  time  integrator.  Diffusive  terms  are  solved  using  a  variant  of  the  implicit  Crank- 
Nicholson  method  (also  second-order  in  time  and  fourth-order  space).  A  corrector  step  is 
then  applied  to  the  predictor  velocities  to  enforce  zero  divergence  at  the  new  time  step.  This 
operation  requires  solving  a  Laplaces  equation,  at  each  time  step.  In  summary,  advancing  a 
solution  one  time  step  requires  solving  four  linear  systems  for  the  fractional  step  velocities 
and  temperature  and  one  Laplaces  equation  for  the  pseudo-pressure  <f>  at  the  next  time  level. 

The  above  algorithm  was  implemented  using  the  PETSc  libraries.  Storage  savings  are 
incorporated  by  only  allocating  storage  for  locally  owned  data  on  each  process  and  the  ghost 
points  from  adjacent  processes.  Similarly,  coefficients  of  linear  systems  are  only  saved  for 
local  grid  points.  A  sparse  storage  technique  has  been  implemented  to  speed  up  computation 
the  solution  of  linear  systems  and  save  storage. 
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The  solution  of  the  fractional  step  velocities  and  temperature  field  converge  quickly  using 
the  PETSc  provided  GMRES  Krylov  subspace  solver.  However,  solving  for  the  pseudo¬ 
pressure  involves  a  much  stiffer  system  due  to  the  full  Neumann  boundary  conditions  nec¬ 
essary  at  the  rigid  boundaries  of  the  domain  for  the  incompressible  N-S  equations.  The 
singularity  of  the  system  is  addressed  by  removing  the  null  space  from  the  solution. 

Two  main  versions  of  the  code  were  created:  a  doubly-periodic  channel  code  with  periodic 
boundary  conditions  in  the  streamwise  and  spanwise  directions,  and  a  streamwise  spatial 
code  with  a  periodic  spanwise  direction.  The  spatial  code  was  modified  for  an  Ekman  layer 
or  channel  through  a  simple  adjustment  of  the  boundary  conditions. 


III.  SUMMARY 

The  objective  of  this  research  was  to  investigate  in  detail  the  dynamics  of  turbulence 
in  simple  and  complex  turbulent  flows,  primarily  with  regards  to  understanding  the  contri¬ 
butions  of  coherent  structures  to  the  turbulence  generation  cycle  and  the  ability  of  closure 
models  to  accurately  approximate  turbulence.  In  ascertaining  the  interactions  and  roles  of 
energetical  structures  as  well  as  their  relationships  with  intercomponent  energy  transfer  and 
overall  turbulent  kinetic  energy,  the  underlying  complex  mechanisms  behind  turbulence  are 
now  better  understood.  The  examination  of  turbulence  and  heat  transfer  closures  also  as¬ 
sisted  this  objective  as  in  addition  to  attaining  significant  reductions  of  computational  costs, 
the  understanding  of  intercomponent  energy  transfer  and  turbulence  production  is  crucial 
within  model  design  and  improvement.  In  parallel,  we  have  applied  passive  flow  control 
by  phononic  subsurfaces  on  three-dimensional  disturbances,  exhibited  by  flows  undergoing 
secondary  (K-type)  transition,  and  demonstrated  that  this  concept  is  successful  even  for 
large- amplitude  nonlinear  instabilities.  The  next  step  in  the  research  is  to  combine  the  two 
thrusts  and  design  phononic-crystal  structures  tuned  for  turbulent  flow  in  order  to  realize 
systemic  suppression  of  targeted  coherent  structures. 

The  publications  referenced  in  Refs.  [5],  [18],  and  [23]  have  been  direclty  funded  by  this 
research. 
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